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We calculate the low temperature properties of an acoustic polaron in three dimensions in thermal 
equilibrium at a given temperature using a specialized path integral Monte Carlo method. In 
particular we find numerical evidence that the chosen Hamiltonian for the acoustic polaron describes 
a phase transition from a localized state to an unlocalized state for the electron as the phonons- 
electron coupling constant decreases. The phase transition manifests itself with a jump discontinuity 
in the potential energy as a function of the coupling constant. In the weak coupling regime the 
electron is in an extended state whereas in the strong coupling regime it is found in a self-trapped 
state. 
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I. INTRODUCTION 

An electron in an ionic crystal polarizes the lattice in its neighborhood. An electron moving with its accompanying 
distortion of the lattice has sometimes been called a "polaron"i£. Since 1933 Landau addresses the possibility whether 
an electron can be self-trapped (ST) in a dcformable lattice 3 - - — . This fundamental problem in solid state physics has 
been intensively studied for an optical polaron in an ionic crystal 6 -! 7 -^— . Bogoliubov approached the polaron strong 
coupling limit with one of his canonical transformations. Feynman used his path integral formalism and a variational 
principle to develop an all coupling approximation for the polaron ground stated. Its extension to finite temperatures 
appeared first by Osak a 13 ' 14 , and more recently by Castrigiano et alJ^—. Recently the polaron problem has gained 
new interest in explaining the properties of the high T c superconductors^. The polaron problem has been also studied 
to describe impurities of lithium atoms in Bose-Einstein ultracold quantum gases condensate of sodium atomsi^. In 
this context evidence for a transition between free and ST polarons is found whereas for the solid state optical polaron 
no ST state has been found yel£~— . The Bogoliubov dispersion at low k is similar to that of acoustic phonons rather 
than to that of optical phonons. 

The acoustic modes of lattice vibration are known to be responsible for the appearance of the ST state d 20 ' 21 . 
Contrary to the optical mode which interacts with the electron through Coulombic force and is dispersionless, the 
acoustic phonons have a linear dispersion coupled to the electron through a short range potential which is believed 
to play a crucial role in forming the ST stated. Acoustic modes have also been widely studied^. Sumi and Toyozawa 
generalized the optical polaron model by including a coupling to the acoustic modes^ 3 -. Using Feynman's variational 
approach, they found that the electron is ST with a very large effective mass and small radius as the acoustic coupling 
exceeds a critical value. Emin and Holstein also reached a similar conclusion within a scaling theory^ in which the 
Gaussian trial wave function is essentially identical to the harmonic trial action used in the Feynman's variational 
approach in the adiabatic limits. 

The ST state distinguishes itself from an extended state (ES) where the polaron has lower mass and a bigger 
radius. A polaronic phase transition separates the two states with a breaking of translational symmetry in the ST 
onei. The variational approach is unable to clearly assess the existence of the phase transition^. Nevertheless Gerlach 
and Loweni concluded that no phase transition exists in a large class of polarons. The three dimensional acoustic 
polaron is not included in this class but Fisher et al£^ argued that its ground state is delocalized. 

In this work we employ for the first time a particular path integral (PI) Monte Carlo (MC) metho d 26 ' 27 to the 
continuous, highly non-local, acoustic polaron problem at low temperature which is valid at all values of the coupling 
strength and solves the problem exactly. Our method differs from previously employed methods 2 ^— since hinges 
on the Levy construction and the multilevel Metropolis method 2 ^. We calculate the potential energy and show that 
like the effective mass it usefully signals the transition between the ES and the ST state. Our results indicate the 
existence of the phase transition. 

The paper is organized as follows: in Sec. [II] we describe the acoustic polaron mathematical model. In Sec. IIIII we 
describe the obscrvables we are interested in. Sec. IIVI contains the description of the numerical scheme used to solve 
the path integral. In Sec. [V]we report our numerical results. And Sec. IVII is for final remarks. 

II. THE MODEL 

The acoustic polaron can be described by the following quasi-continuous model 7 -^, 

A = £ + E ^k h k + E ^ k a k e* k * + H.c.) . 
k k 

Here x and p are the electron coordinate and momentum operators respectively and is the annihilation operator 
of the acoustic phonon with wave vector fc. The electron coordinate a; is a continuous variable, while the phonons 
wave vector k is restricted by the Debye cut-off k . The first term in the Hamiltonian is the kinetic energy of the 
electron, the second term the energy of the phonons and the third term the coupling energy between the electron and 
the phonons with an interaction vertex Tk = huk {S / po) 1 / 2 {k / ko) 1 / 2 where S is the coupling constant between the 
electron and the phonons and po the number density of unit cells in the crystal (po = (47r/3)(fc /27r) 3 in the Debye 
approximation). The acoustic phonons have a dispersion relation LUk = uk, u being the sound velocity. 

Using the path integral representation (see ReL^ 2 - section 8.3), the phonon part in the Hamiltonian can be exactly 
integrated owing to its quadratic form in phonon coordinates, and one can write the partition function for a polaron 
in thermal equilibrium at an absolute temperature T (/? = l/fe^T, with ks Boltzmann constant) as follows, 

dx / / e-T, s ^)^^Vx(t), (1) 

J-Jx=x(o) 
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eff 



where the action S is given b y 35 ; 36 

s = t r - th r *r * / «<>-*<•»-"-' 

= S f +U. (3) 

Here 5/ is the /ree particle action, and W the inter-action and we denoted with a dot a time derivative as usual. 
Setting h = m ~ uk = kg = 1 the inter-action becomes, 

W= / dt dsV eff (\x(t)-x(s)\,\t-s\), (4) 
Jo Jo 

with the electron moving subject to an effective retarded potential. 

-— dggeV^ 1 (5) 

where q = k/k , Id = / ? < 1 ^9 = 47r/3, and we have introduced a non-adiabatic parameter 7 defined as the ratio of 
the average phonon energy, huk , to the electron band-width, (hk ) 2 /2m. This parameter is of order of 10 -2 in typical 
ionic crystals with broad band (~eV) so that the ST state is well-defined^ 3 -. In our simulation we took 7 = 0.02. 
One can expect two kinds of polarons: Electrons in alkali halides and silver halides are nearly ES while holes in alkali 
halidcs are in the ST stated. The hole is ES in AgBr— and ST in AgCl 3 ^. The most dramatic observation of the 
abrupt change of exciton from ES to ST states was made on mixed crystals AgBri-zCLj—S. 

III. THE OBSERVABLES 

The free energy (the ground state energy, E, in the large j3 limit) of the polaron is F = —{dZ/df3)/Z = (JC + V), 
where the first term is the kinetic energy contribution, /C, and the second term is the potential energy contribution, 
V . We have, 

Scaling the Euclidean time t = (3t' and s = (3s' in Eq. @, deriving U with respect to f3 and undoing the scaling in 
the end we obtain for the potential, 



SS ^ f p f 1 sm(J + q\x(t)-x(s)\ 

V = -^R dt ds d 11^lf- J -e-^(2~q\t- S \). 

2 P Jo Jo Jo Jlq\x(t)-X{s)\ 



Taking the derivative with respect to (3 of the action after having scaled both the time as before and the coordinate 
x = and undoing the scaling in the end we obtain for the kinetic energy, 



3S f f f 1 
K, = —— I dt ds dqq 



sin U /|g|a:(i) -x(s)|) 

cos I J-q\x(t)-x(s)\ - 



4/3 7o 

In the following we will be concerned with a numerical determination of the potential energy. 



Jlq\x(t)-x(s)\ 



e -g\t-s\ 



(7) 



IV. PATH INTEGRAL MONTE CARLO 

To calculate the PIMC, we first choose a subset of all paths. To do this, we divide the independent variable, 
Euclidean time, into steps of width r = (3/M. This gives us a set of times, tk = kr spaced a distance r apart between 
and (3 with k = 0, 1, 2, . . . , M. At each time tk we select the special point Xk = x(tk), the /cth time slice. Wc 
construct a path by connecting all points so selected by straight lines. It is possible to define a sum over all paths 
constructed in this manner by taking a multiple integral over all values of x^ for k = 1, 2, . . . , M — 1 where x = x a 
and X]\j = xi, are the two fixed ends. The simplest discretized expression for the action can then be written as follows, 

M , , 2 M M 

S = £ {Xk -'~ Xk) + r 2 E E V(U, (8) 

k=l i=l j = l 
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where V(U,tj) = V ef f(\xi 



■ j\) is a symmetric two variables function. In our simulation we tabulated this 



function taking |aj^ — Xj\ = 0, 0.1, 0.2, . . . , 10 and \i — j\ =0,1,..., M. The total configuration space to be integrated 
over is made of elements s = {xq,Xi, . . . where Xk are the path time slices subject to the periodic boundary 

condition xm = xq. In order to compute the potential energy P = (V) in the simulation we wish to sample these 
elements from the probability distribution, 7r(s) = e /Z, where the partition function Z normalizes the function ir 
in this space. 

In our simulation we chose to use the bisection method, a particular multilevel MC metho d 26 ! 41 ' 42 , with correlated 
sampling. The transition probability for the first level is chosen as T\ oc exp[(xj +TO /2 — x) 2 /2a 2 (m/2)} where m = 2', 
I being the number of levels, x = (x, + Xi +m )/2 and a(to/r) = <J ([x(t) — (x(t + to) + x(t — to))/2] 2 ) (for the 
first levels these deviations are smaller than the free particle standard deviations used in the Levy construction 4 ^ 
o~f(£) = \/£t/2 with £k = m/2 h in the fcth level. Much smaller in the ST state.). And so on for the other levels: 
s 2 = {x l+m / 4 ,x l+3m / 4 }, ... , si = {x i+ i,x i+2 , ■ ■ ■ ,x i+m _i}. And s = {aJo, . . . , x it x i+m , . . . , xm-i} where i is 
chosen randomly. Calling tt(s) — eT u jZ ', the level inter-action is 7Tfc(soj ■ • ■ , Sfc) = / dsk+i ■ ■ ■ dsi tt(s). For the fcth 
level inter-action we thus chose the following expression, 



7Tfc cx exp 



[M/i k ] [M/e k ] 
-(^4) 2 J] E V(i£ k r,j£ k r) 



i=l 



(9) 



In the last level 



1 and the level inter-action 7n reduces to the exact inter-action jr. The acceptance probability for 



the first level will then be, A\ = min 



P 1 (s) #i(s')7r (s) 
' P 1 (s') #i(s)#o(s') 



with Pi oc exp{— (x i+m / 2 — x) 2 {l/o- 2 (m/2)- 



l/a 2 (m/2)]/2}. 



The initial path was chosen with all time slices set to 0. During the simulation we maintain the acceptance ratios 
in [0.15,0.65] by decreasing (or increasing) the number of levels in the multilevel algorithm as the acceptance ratios 
becomes too low (or too high). We will call Monte Carlo step (MCS) an attempted move. 



V. RESULTS 



We simulated the acoustic polaron fixing the adiabatic coupling constant 7 = 0.02 and the inverse temperature 
,3=15. Such temperature is found to be well suited to extract close to ground state properties of the polaron 3 ^. For 
a given coupling constant S we computed the potential energy P extrapolating (with a linear \ square fit) to the 
continuum time limit, r — > 0, three points corresponding to time-steps choosen in the interval r <G [1/100, 1/30]. In 
Fig. [T]and Tab. [J we show the results for the potential energy as a function of the coupling strength. It is clear the 



TABLE I. MC results for P as a function of S at /3 = 15 and 7 = 0.02 displayed in Fig. □ The runs where made of 5 x 10 s 
MCS (with 5 x 10 4 MCS for the equilibration) for the ES states and 5 x 10 6 MCS (with 5 x 10 5 MCS for the equilibration) for 
the ST states. 
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10 


-0.573(8) 


20 


-1.17(2) 


30 


-1.804(3) 


40 


-2.53(3) 


50 


-3.31(4) 


53.5 


-3.61(1) 


55 


-11.4(3) 


60 


-16.1(5) 


70 


-23.3(3) 


80 


-30.0(3) 



transition between two different regimes which correspond to the so called ES and ST states for the weak and strong 
coupling region respectively. We found that paths related to ES and ST are characteristically distinguishable. Two 
typical paths for the ES and ST regimes involved in Fig. Q] are illustrated in Fig. [2] The path in ES state changes 
smoothly on a large time scale, whereas the path in ST state do so abruptly on a small time scale with a much 
smaller amplitude which is an indication that the polaron hardly moves. The local fluctuations of the x(t) and of the 
potential energy 'P(MCS) have an auto-correlation function which decay much more slowly in the ES state than in 
the ST one. Moreover the ES simulations are more time consuming than the ST ones. 
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Concerning the critical property of the transition between the ES and ST states our numerical results are in favor 
of the presence of a discontinuity in the potential energy. Even if there is no trace of a translational symmetry 
breaking as shown by the ST path in Fig. [2] where the initial path was x{t) = for all t. With the increase of 
/3, the values for the potential energy P increase in the weak coupling regime but decrease in the strong coupling 
region. From second order perturbation theory (see RefJ^ section 8.2) follows that the energy shift 2?(7, 5) is 
given by — 357[l/2 — 7 + 7 2 ln(l + 1/7)] from which one extracts the potential energy shift by taking P(7, 5) = 
7<iE(7, S)/dj. From the Feynman variational approach of Rcf^ follows that in the weak regime the energy shift is 
— 357[l/2 — 7 + 7ln(l + I/7)] and in the strong coupling regime 
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FIG. 1. Shows the behavior of the potential energy P as a function of the coupling constant S. The points are the MC results 
(see Tab. U, the dashed line is the second order perturbation theory result (perturbation) and the dot-dashed line is the 
variational approach from Ref.— (variational) in the weak and strong coupling regimes. 
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FIG. 2. The polaron closed path x(t) = (X (t) ,Y (t) , Z(t)) as a function of Euclidean time t in units of r. 



Note that since S and r appear in the combination St 2 in U (and St in V) the same phase transition from an ES 
to a ST state will be observed increasing the temperature. With the same Hamiltonian we are able to describe two 
very different behaviors of the acoustic polaron as the temperature changes. 



VI. CONCLUSIONS 



In conclusion we used for the first time a specialized PIMC method to study the low temperature behavior of an 
acoustic polaron. At an inverse temperature j3 = 15, close to the ground state of the polaron, and at a non-adiabatic 
parameter 7 = 0.02, typical of ionic crystals, we found numerical evidence for a phase transition between an ES state 
in the weak coupling regime and a ST one in the strong coupling regime at a value of the phonons- electron coupling 
constant S* = 54.3(7) in good agreement with the prediction of Refi^, S* ~ 1/7, and the MC simulations of Refi22. 

To understand the motion of an electron in a deformable lattice we have to consider the fact that the interaction 
of the electron with the acoustic phonons induces a well barrier proportional to the coupling constant, but it decays 
due to the retarded property. In the weak coupling region, the electron can easily tunnel through the barrier so that 
it almost freely moves in the lattice. One can regard the tunneling of the electron as an indication of the motion of 
the polaron. In this case, a few phonons arc involved and the acoustic polaron has a small mass which has a similar 
magnitude to the mass of the electron and a large radius. In the strong coupling region, the well barrier becomes 
sufficiently deep and the electron is temporarily bounded in the polaron and cannot tunnel through the barrier until 
it gains enough energy from the phonons. Much more phonons are involved in the polaron and the mass of the 
polaron becomes much greater than that of the electron with a small radius. In this argument the specific form of 
the interaction vertex is of fundamental importance. 

We used for the first time a PIMC with the bisection method and correlated sampling as an unbiased numerical 
mean to probe the low temperature properties of the acoustic polaron. This is an independent route to the Fcynman's 
variational approach and proved to give reliable results on the existence of the ST state for the electron in a deformable 
lattice as conceived by Landau. However, the self-trapping we observe in our numerical analysis is not a complete 
localization of the electron within the polarization cloud of the phonons. The electron path still undergoes small, 
nearly uncorrelated, fluctuations in Euclidean time in this ST state. Our numerical results support the presence of 
a discontinuity in the potential energy as a function of the coupling constant and this would be an indication of the 
existence of a phase transition between the ES and the ST states even if we found no trace of a translational symmetry 
breaking. Moreover The discontinuity of the ES/ST transition, if it exists, may depend on the cutoff parameter as 
pointed out by e.g. Refi^. In the cold atoms context the role of this parameter is important and the existence of 
a discontinuous transition is more questionable than in the solid state polaron case^. The present study reports a 
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single value study of the acoustic polaron case thus restricting the generality of the conclusions. 

In a truly localized state the polaron should not diffuse at all (strong localization) or at least should attain a 
subdiffusive behavior (weak localization). But these properties can be checked by looking at the real time dynamics 
of the system and cannot be checked by the Monte Carlo methods as that used in this work which deals with polaron 
properties in imaginary time. Clearly our numerical MC results support the claim of a dicontinuity but do not give 
any proof of a localization transition. 

A possible further study could involve the dynamic properties associated with the two different types of motion 
and bipolarons for short range interacting systems. 
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